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Mountain glaciers respond directly to changes in precipitation and temperature, thus their margin extent is 
a high-sensitivity climate proxy. Here, we present a robust ^*^Be chronology for the glacier maximum areal 
extent of central Spain paleoglaciers dated at 26. 1 ± 1 .3 ka BP. These glaciers reached their maximum extent 
several thousand years earlier than those from central Europe due to the increased precipitation within a 
cold period between 25 to 29 ka BP, as confirmed by a local speleothem record. These paleoclimate 
conditions impacted the maximum extent of mountain glaciers along the western and central 
Mediterranean region. The cause and timing of the enhanced precipitation implies a southward shift of the 
North Atlantic Polar Front followed by storm tracks in response to changes in insolation via orbital 
parameters modulation. Thus, these mountain paleoglaciers from the Mediterranean region record an 
ocean-continent climate interaction triggered by external forcing. 



The Mediterranean region has a westerly circulation due to its mid-latitude position, and therefore, its climate 
is directly impacted by the North Atlantic. During the last glaciation, the mid-latitudes of the North Atlantic 
witnessed the most significant changes in the surface ocean dynamics \ which greatly affected the global 
climate. Therefore, the Mediterranean mountain paleoglaciers^'^ record the ocean- continent climate interactions 
during the last glaciation in a critical location. However, there are only a limited number of locations with robust 
chronologies for the glacier maximum extent (GME) in the Mediterranean mountains to allow the identification 
of regional paleoclimate patterns (Fig. 1, Supplementary Table SI). 

The paleoglaciers from the Spanish Central System have an ideal geologic and geomorphologic setting to 
provide a robust chronology for the regional GME (see supplementary information). Additionally, a U-Th dated 
paleoclimate record based on speleothems from the nearby Eagle Cave, provides information on the amount of 
precipitation supplied to the region during the period of interest, allowing our study to evaluate the impact of 
precipitation on the areal extent of the glaciers. 

Many summits and plateaus of the Spanish Central System hosted glaciers during the last glacial period (Fig. 1), 
although they are absent today. The mountains are mostly composed of granite and gneiss rocks, which provided 
tills and erratics with boulders of great size^ ideal for using the cosmogenic ^°Be exposure dating technique. The 
paleoglaciers from the region have a morphosedimentary record with differentiated homogeneous geomorpho- 
logic units^ The "Peripheral Deposits" (PD) is the unit that contains the GME, although also records later pulses 
or stages. The PD consists of a series of geomorphic indicators (i.e., moraines, boulders-belts and erratic strewn 
boulders). Frequently, this unit has a sequence of indicators (e.g., several parallel moraines) whereas the 
most external indicators represent the GME and the internal indicators record later fluctuations of the ice 
margin in which the glaciers were close to the GME extension (see supplementary information for further 
descriptions and illustrations of the PD). The overflow of ice on top of the plateaus or gentle slopes aside from 
fluvial or slope dynamics has favored the complete preservation of PD sequences in many lateral moraines. 
Therefore, recognition of these sites after detailed geomorphologic mapping minimizes the risk of misidentifica- 
tion of the GME. 
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Figure 1 | Map of the Spanish Central System showing the limits of the paleoglaciers of the region during their maximum extension. Detailed sketches 
show the sampling sites for ^°Be. The inset map of the Mediterranean mountains shows the location of different sites discussed in the text: 1) MD95- 
2040'^ 2) MD95-2042^^ 3) Villars Cave'\ 4) Taghamento amphitheater^^ 5) central Apennines'^ 6) Uludag Mnt.'', 7) Sandiras Mnt/°, 8) Erciyes Mnt.". 
Mapping of Bejar and Credos glaciers after Pedraza et aV. The digital elevation model used in this map is from the free access database of the ICN 
(www.ign.es). 



Results 

Nine paleoglaciers from the Bejar, Credos and Guadarrama ranges 
were studied to determine the GME chronology and its replicability. 
Twenty-five samples were obtained from boulders deposited by the 
paleoglaciers in the PD unit (Supplementary Table S2). To minimize 
the bias of the age determination and to improve the accuracy of the 
^°Be chronology, a specific erosion rate was applied to the age cal- 
culation for each sample. However, the effect of erosion on the cal- 
culated dates assuming no erosion differs less than 1.0% in relation to 
the reported dates using the modeled erosion rate (Supplementary 
Table S2). The samples were always collected from geomorphic indi- 
cators representing glacier edges. So, the dates reported are inter- 
preted as ages of sedimentation that records the timing of the glacier 
pulse that reached that particular location. The dating results show 
two age clusters for the PD unit, which are confirmed to be statist- 
ically different subpopulations according to Kolmogorov-Smirnov 
test (p = 0.00041). Therefore, due to the clear statistical separation 
between the two age subpopulations and the occurrence of major 
climate events between both periods, we recognize two different 
glacier stages (Fig. 2; Fig. 3). The age of each stage results from the 
highest probability after the sum of the probability distribution func- 
tions (PDFs) of every sample, and the reported uncertainty is the 
standard error of the mean of the two compiled PDFs. The stage 
representing the GME has a depositional age of 26.1 ± 1.3 ka BP, 
and is found in the most external geomorphic indicators. The second 
stage has an age of 2 1 .3 ± 0.7 ka BP and contains boulders from both, 
internal and external geomorphic indicators (Fig. 2; Fig. S2-S4). The 



chronology of the PD unit confirms the existence of an early GME in 
central Spain (n = 6 in 4 paleoglaciers) in comparison with glacier 
chronologies from central Europe: the glaciers from the northern 
Alps^, the southern sector of the Scandinavian Ice Sheet^ or the 
Great Britain-Irish Ice Sheet^. However, the second glacier stage 
recorded in central Spain (n = 16 in 8 glaciers) is in agreement with 
the central Europe and eastern Mediterranean GME glacier chro- 
nologies^"^ ^ and with secondary glacier stages in the mountain gla- 
ciers from the southern Alps^^ and the Apennines^^'^^. The PD unit 
has a limited extent (averaging <250 m along sequence traverses), 
and most of the record was accumulated during the later stage, which 
implies that during both stages the extents of the glaciers were sim- 
ilar. Thus, some geomorphic indicators identified as external have 
ages that do not correspond with the central Spain GME, which 
implies that boulders from a further glacier extent have not been 
preserved/identified on these paleoglaciers or that the ice margin 
of the later stage overrode the previous deposits at that particular 
location, illustrating the limited extent differences in some cases 
between both stages. 

Two stalagmites from Eagle Cave, located only 10 km south of the 
Credos Range, were U-Th dated (Supplementary Table S3) and their 
isotope signal was studied. The 5^^0 record of the samples depends 
on several factors and its variability cannot be interpreted as the 
result of a single control (see supplementary information). 
However, most of the 5^^0 variability depends on the amount of 
precipitation with the changes in moisture source (mostly account- 
ing for differences in the seasonal ice cover in the North Atlantic), 
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Figure 2 | Probability distribution functions (PDFs) of the ^"Be dates in PD unit of paleoglaciers from Bejar, Credos and Guadarrama ranges, (a) PDFs 
of the ^"^Be dates. Samples from each paleoglacier have different colours. Dashed lines represent discarded dates, (b) Cumulative PDFs for boulders 
(bold lines) from internal and external geomorphic indicators. Note that the dates of the external geomorphic indicators provide two subpopulations 
(Ext-A and Ext-B). PDFs of individual samples are also plotted, (c) Similar to (b) but divided in three panels in order to differentiate the individual 
samples in each subpopulation. The black vertical lines are the most probable age for each subpopulation and the grey bars are the standard error of the 
mean for each cumulative PDF. (d) PDFs of the two stages recorded in the PD unit (bold lines) together with the individual PDFs. 



likely playing a secondary but significant role. Additionally, the min- 
eralogy is a factor affecting the 5^^0 signal during two particular 
periods around 24 and 30 ka BP when aragonite, instead of calcite, 
comprises most of the carbonate (Fig. 85). Thus, lower values in the 
5^^0 record are interpreted as wetter periods, with moisture source 
potentially contributing in part to the variability of this proxy (Fig. 3). 
During the Heinrich events H2 and H3, the climate in the region was 
extremely dry^^ as reflected by high 5^^0 values and the presence of 
aragonite. The lower values of the 5^^0 record between 25 to 29 ka 
BP when compared to the period 19 to 23 ka BP, in which the central 
Spain paleoglaciers record the second glacier stage, confirm that the 
climate was wetter in this region at the time of GME than during the 
second glacier stage. A similar paleoclimate pattern is recorded in 
other sites of the Iberian and the Italian peninsulas^^'^^ supporting its 
regional influence over part of the Mediterranean. 

Discussion 

The extent of mid-latitude mountain glaciers depends on their long 
term mass balance. Therefore, glacier extent changes are related to 
persistent variations in their accumulation and/or ablation values. In 
general, the accumulation is primarily controlled by the cold season 
snow precipitation, while the ablation depends mostly on the warm 
season temperature, controlled in part by the energy received by the 
sun and modulated by the orbital parameters. At 40°N latitude, the 
integrated summer insolation (ISI) that affects the glacier ablation 



had a minimum at 27 ± 4 ka^^. However, sea surface temperature 
(SST) offshore Iberia in the Atlantic Ocean was cooler during the 
period 19 to 23 ka BP than during the central Spain GME^^. Even 
lower SSTs were recorded during Heinrich events, but particularly 
dry conditions associated with these periods prevented the glaciers 
in southern Europe from reaching their GME^. In fact, the cause for 
an early GME in the Mediterranean region in comparison central 
Europe is thought to be related to a glacier mass balance dominated 
by increased precipitation rather than cooler temperatures^^. 
Nevertheless, temperatures in southern Europe were so cold by 
30 ka that in Villars Cave (southern France), the speleothems stop 
their growth due to continuous permafrost over the cave^\ This local 
and regional paleoclimate context suggests that the optimal com- 
bination of wet and cold conditions that caused the paleoglaciers 
to reach their maximum extension were recorded in central Spain 
around 25-29 ka BP. Thus, the cooler and drier climate conditions 
recorded between 19 and 23 ka BP in comparison with the 25-29 ka 
BP period were responsible for glaciers of large dimensions in the 
region, but not enough for reaching the GME. 

In order to investigate the causes of the wet period between 
Heinrich events H2 and H3 and to discard the existence of other 
periods suitable to record the GME, the regional paleoclimate con- 
text is shown in figure 4 up to the onset of the last glaciation at the end 
of marine isotope stage 5 (MIS5). During the 19-23 ka BP period, the 
North Atlantic Polar Front (PF) was at the latitude of Iberia^^. The 
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Figure 3 | Chronology of the PD in the Central Spain paleoglaciers and 
paleoclimate setting, (a) Probability distribution functions of the ^°Be ages 
for the two stages. Black vertical lines represent the stage ages with the grey 
vertical bars being the standard error of the mean, (b) ISI variability for 
40°N^^ for most likely x configurations; x = 325 (dashed); 350 (solid) and 
375 (dotted), (c) Eagle Cave 5^^0 record (light blue) and ice volume 
corrected 5^^0 record (dark blue) with a 0.2 ka filter (bold blue). The 
isotope ratio correction takes into account sea level changes by Peltier and 
Fairbanks^^ considering a rate of change of 0.06%o/10 m of sea leveP^. (d) 
Black rectangles labelled HI to H3 refer to Heinrich events^*^. (e) Green 
rectangle corresponds with the period of lack of calcite deposition in Villars 
Cave (southern France) due to extreme cold conditions^ \ Error bars 
indicate U-Th ages uncertainty. 

SST difference of the Ocean cores MD95-2040 and MD95-2042 
located along the coast of Iberian margin in the North Atlantic 
reflects the proximity of the PF, where the thermal gradient is max- 
imum\ During the past 80 ka, two periods of outstanding southern 
displacement of the PF are recorded: between 26-32 ka BP and 66- 
72 ka BP. The storm tracks over Europe followed these southern 
shifts since their energy is taken from the ocean thermal gradient^^'^^. 
The resulting climate setting provides more precipitation over most 
of the Mediterranean region, with the exception of the periods of 
extreme low SST during Heinrich events, which limited the moisture 
uptake^^. Both periods of increased thermal gradient along the 
Iberian coast roughly match the minima in the ISI at 40 °N, and 
are candidates to record the GME in the Mediterranean region. 
However, the continental temperatures in the early part of MIS4, 
when ISI conditions are minimal (72 ± 4 ka BP) were too mild 
compared with the later period (27 ± 4 ka BP), when conditions 
were then optimal to reach the GME in the region. 

While the glacier records in central Europe reach their GME 
within the 19-23 ka period^"^, paleoglaciers from central Iberia, 
southern Alps^^ and central Italy^^'^^ record an early GME around 
26 ka BP. In contrast, paleoglaciers from Turkey have shown a GME 



during the 19-23 ka BP period^"^\ This difference is attributed to the 
intensification of the Siberian High during glacial times^^, which 
prevented the southward shifted storm track system to significantly 
impact continental Asia in the eastern Mediterranean, maintaining 
similar precipitation conditions between 20 and 30 ka BP in the 
region^^. 

In the western and central Mediterranean, a persistent increase in 
the snow accumulation during the cold period from 25 to 29 ka BP 
caused mountain glaciers to reach their GME. This early GME in 
relation to glaciers from central Europe and the eastern Mediter- 
ranean is the result of differential precipitation distribution in the 
continent. The increased precipitation over part of the Mediter- 
ranean was the consequence of a southward shift of the storm tracks 
in response to ocean dynamics in the North Atlantic. The match of 
the limited seasonal amplitude determined by the lower obliquity 
and normalized orbital forcing^^'^^ with the larger SST gradients 
along Iberian coast, suggests a substantial role of the winter sea ice 
cover in controlling the precipitation regime over the Mediterranean 
region. The paleoglaciers from central Europe and the eastern 
Mediterranean record their GME between 19 and 23 ka BP, a cold 
and relatively dry period in which mountain glaciers from the west- 
ern and central Mediterranean also record important glacier stages, 
although these do not reach their GME. In any case, in all these 
regions the GME took place during the MIS2, fitting the period of 
minimum atmospheric CO2, temperatures and relative sea level that 
commenced shortly after 30 ka BP^^. 

Methods 

Estimation of the erosion rate of boulder samples. We conducted a study of the 
regional erosion rates in mountains of central Spain in order to minimize the 
potential bias of the calculated ^°Be exposure ages in boulders deposited by glaciers 
and to improve accuracy of the chronology. To further constrain the erosion since 
deglaciation, we carried out calculations based on geomorphic evidence. Preservation 
of striations and polished surfaces due to glacial activity are uncommon in mountains 
from central Spain, but grooves and other signs of glacier abrasion are well 
represented. In these glacially abraded surfaces, the feldspar phenocrysts and quartz 
veins are resistant features due to differential weathering. Faint striations on top of 
these outstanding features are preserved under optimal conditions, supporting the 
assumption that their weathering since deglaciation is negligible. Therefore, the 
outstanding height of these features was considered to measure the erosion since 
deglaciation. Erosion was measured in ten sites along Bejar, Credos and Guadarrama 
ranges. Selected sites are in top sections of roches moutonnees where concentrated 
erosion due to water runoff or fluvial activity is unlikely. Ponds or other sites of 
potential rainwater accumulation were avoided. Average erosion in each station was 
considered as the mean of fifty measures in different spots. To calculate the erosion 
rate a deglaciation age for each sampling site was considered. The deglaciation stage of 
the studied sites was established after regional mapping^ and available regional 
chronologies^''"^^. Two different deglaciation stages with arbitrary deglaciation ages of 
15 ± 1 and 19 ± 1 ka cover the location of all the studied sites. Although age 
determination for each stage introduces a significant uncertainty, the constraint of 
deglaciation on these mountains to short periods makes erosion rates less vulnerable 
to age uncertainties. One of the sites shows high erosion rate (3.58 ± 1.27 mm*ka~^) 
in agreement with field observations. The abundant presence of xenoliths in these 
rock suggest that its particular petrologic properties make it prone to heavy 
weathering in comparison with the rest of granites and gneisses analyzed. This 
observation is corroborated by tests done with a Schmidt hammer that provided Q- 
values of 37.4 ± 2.8 (dimensionless) on the heavily weathered boulders in comparison 
with the rest of analyzed samples with Q-values from 47.6 ± 2.1 to 53.2 ± 2.1. These 
heavily weathered boulders are easily recognizable in the field and they were 
discarded for ^°Be exposure dating sampling. For the rest of the sites studied, the 
erosion rate was between 0.58 ± 0.16 and 0.29 ± 0.08 mm*ka"\ These low erosion 
rates for the studied mountain ranges are in good agreement with previous 
weathering studies in weathering pits^^ and the preservation of abrasion features since 
deglaciation. There is a positive correlation of erosion rate and elevation (r^ = 0.94; p 
< 0.01) following a linear relationship (Fig. S6). Therefore, this relationship was used 
to model the erosion rate for every sampled boulder and correct the ^°Be exposure 
ages. Although this method for estimating erosion rates has inherent limitations, it 
improves the accuracy of the ^°Be exposure ages in comparison with traditional 
assumptions for erosion rate. Thus, we avoid assuming an unrealistic zero erosion 
since deglaciation, or generalizing regional erosion rates. 

Calculation of ^"Be dates and chronology of PD stages. The samples collected from 
PD unit followed standard procedures of sample preparation for ^"Be dating^"^. The 
^°Be/''Be ratio was measured at PRIME laboratory. Sample shielding factors and 
sample thickness were measured in the field and incorporated in the calculations. We 
suspect that boulders remained windswept due to their height, and so snow cover 
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Figure 4 | Paleoclimate context affecting paleoglaciers in the Mediterranean region, (a) SST difference between cores MD95-2042 and MD95-2040^'^ at 
2 ka intervals (thin red) with a 6 ka fdter (bold red), (b) SST at MD95-2042^^ (c) Obhquity"^ (d) Normalized orbital forcing"", (e, f) ISI variability and 
period of calcite deposition in Villars Cave as in figure 3. Heinrich events HI to H6 are represented as vertical grey bars^'^. 



estimates were not considered. The coordinates of samples were taken with a GPS in 
the field, and their location was improved by relocating the boulders in a 5 m cell 
resolution digital elevation model with 0.25 m pixel size resolution orthophoto 
incorporated, in which individual boulders were identified. Calculation of the 
production rate and ages was done with the CRONUS-Earth online calculator v2.2^^ 
Because there are no production rates for the specific area of our research, we have 
used the standard scaling and accepted ^°Be production rates in Balco et aV^ and 
chosen to report the external uncertainty given the small discrepancies between the 
scaling methods available. Our results and conclusions regarding the early GME and 
two phases of glacial activity are insensitive to which of the scaling scheme we choose. 
Here, we report our data based on Lifton et aV^ scaling scheme given in the 
CRONUS-Earth online calculator that includes temporal variations in the cosmic ray 
flux. Ages are provided in ka BP (before present), where present is considered the year 
1950 (i.e., ages provided by the online calculator were shifted 58 to 61 years depending 
on the sampling year) in order to match other chronologies. Two sets of ages are 
displayed in Table S2: considering no erosion rate and the modelled erosion rate. The 
chronology of this paper consider the modelled ages, although differences are < 1% in 
relation to those that consider zero erosion rates. 

The modelled age results and the external uncertainties (la) were used to calculate 
the probability distribution functions (PDFs) of each date. From the 25 dates calcu- 
lated 2 of them (PE-25 and HG-65) were considered unreliable due to lack of over- 
lapping (within la external uncertainty) with samples from the same geomorphic 
indicator. Most likely, these samples suffered some unexpected rotation or surficial 
erosion after exposition that made their ages younger than other boulders in the same 
paleoglacier sequence. These dates are discarded in further calculations. 

The distribution of dates of PD unit shows two clear subpopulations in the PDFs 
(Fig. 2A). The dates are categorized based on the geomorphologic context of the 
sampled boulders (internal or external geomorphic indicators) along the sequences of 
the PD unit, and the cumulative PDFs have been calculated using k- means clustering 
to split the data (Fig. 2B and Fig. 2C). The internal geomorphic indicators provided 
dates within one subpopulation with most probable age at 20.1 ka BP ("Int" in Fig. 2). 
The external geomorphic indicators provided dates distributed in two subpopula- 
tions, with most probable ages at 26.1 and 22.2 ka BP ("Ext- A" and "Ext-B" 
respectively in Fig. 2). The Kolmogorov-Smirnov test was applied to evaluate if these 
populations are statistically different. Thus, the two subpopulations found on the 
external geomorphic indicators (Ext-A and Ext-B) belong to two different popula- 
tions (p = 0.00036), whereas the younger subpopulation of dates from the external 
geomorphic indicators (Ext-B) and the dates from the internal geomorphic indicators 
(Int) can be considered to belong to the same population (p = 0.1441). Therefore, the 
dates from the younger subpopulation of external geomorphic indicators (Ext-B) and 
the dates from internal geomorphic indicators (Int) were merged in a single 



population with most probable age of 21.3 ka BP (Fig. 2D), which is statistically 
different from the subpopulation of most probable age of 26.1 ka BP (p = 0.00041). 
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